function check_limits(CONSTS, plot_data)

    c     = CONSTS.c;
    k0    = CONSTS.k0;
    a     = CONSTS.a;
    d     = CONSTS.d;
    use_V = CONSTS.use_V;

    zeta_step = 0.05; % 0.05
    zeta_max = 1;
    zeta_vec = (zeta_step : zeta_step : zeta_max)'*d;
    int_limit = 10.0;
    integ1_vec = zeros(size(zeta_vec)); % q from int_limit to infty
    integ2_vec = zeros(size(zeta_vec)); % q form 0 to int_limit
    %integ_ch_min_vec = zeros(size(zeta_vec));
    %integ_ch_max_vec = zeros(size(zeta_vec));
    for i = 1:size(zeta_vec, 1)
        zeta = zeta_vec(i);
        m = 1; % 0, 1, 2, 3, 4, ..., 20, ...
        period = pi/(k0*a);
        dots_per_period = 100; % 200
        order = 3;
        q_high = 0.5*log(10)*order/(k0*zeta);
        q1_vec = (int_limit : period/dots_per_period : q_high)';
        q2_vec = (0.0 : period/dots_per_period : int_limit)';

        y1_vec_1 = func_1(q1_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        y1_vec_2 = func_2(q1_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        integ1_1 = trapz(q1_vec, y1_vec_1);
        %integ1_1_g = quad_gauss(@(q)func_1(q, zeta, m, CONSTS), ...
        %    int_limit, q_high, period/4, 4);
        integ1_2 = trapz(q1_vec, y1_vec_2);

        y2_vec_1 = func_1(q2_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        y2_vec_2 = func_2(q2_vec, zeta, m, CONSTS); %zeros(size(q_vec));
        integ2_1 = trapz(q2_vec, y2_vec_1);
        integ2_2 = trapz(q2_vec, y2_vec_2);    

        integ1_vec(i) = - ((2*pi*k0^2*a)/c) * integ1_1 ...
                       + (2*pi/(c*a)) * integ1_2;

        integ2_vec(i) = - ((2*pi*k0^2*a)/c) * integ2_1 ...
                       + (2*pi/(c*a)) * integ2_2;           

        %integ_ch_min_vec(i) = (2*k0/c) * integ_check_for_K1_1(zeta, m, CONSTS) + ...
        %                      (2/(c*k0*a^2)) * integ_check_for_K1_min_V(zeta, m, CONSTS);
        %integ_ch_max_vec(i) = (2*k0/c) * integ_check_for_K1_1(zeta, m, CONSTS) + ...
        %                      (2/(c*k0*a^2)) * integ_check_for_K1_max_V(zeta, m, CONSTS);

    end   

    if (plot_data)
         figure; plot( ...
                      zeta_vec/d, integ1_vec, 'r.-', ...
                      zeta_vec/d, integ2_vec, 'b.-' ...
                      );
         title('K_m(\zeta)');
         xlabel('\zeta /d'); ylabel('K_m'); %ylim([0 8e-13]);
         legend( ...
             ['Numerical integral (q from ', sprintf('%g', int_limit), ' to \infty)'], ...
             ['Numerical integral (q from 0 to ', sprintf('%g', int_limit), ')'], ...
             'location', 'SouthEast');
    end

end